#!/bin/bash
set -e

function info() {
echo Usage: `basename $0` '-m(somatic) -g(single_cell) -l(interval.bed) snp.vcf indel.vcf'
exit 65
}

while getopts ":w:p:l:mgh" opts
do
case $opts in
	w) win_size=$OPTARG;;
	m) somatic=T;;
	g) single_cell=T;;
	l) interval=$OPTARG;;
	p) out_prefix=$OPTARG;;
	\?) info;;
esac
done
shift $(($OPTIND - 1))

test $# -lt 2 && info


. /mnt/ilustre/app/medical/tools/.var


vcf_path=${data_path}/vcf/gatk/
echo
echo $interval
echo dbsnp version is: $dbsnp_version
echo $snpeff_db_version

	
format_af.sh -p$out_prefix.1 $1
format_af.sh -p$out_prefix.2 $2 
anno_info_vcf.pl $1 $out_prefix.1.af.txt > $out_prefix.snp.af.vcf
anno_info_vcf.pl $2 $out_prefix.2.af.txt > $out_prefix.indel.af.vcf


IFS=$'\n'
for dat in $(less $tools_path/script/filters_snp.txt); do

filter=$(echo $dat |cut -f2)
filter_str=$(echo $dat|cut -f1)
filter0=$(echo --filterExpression \"$filter\" --filterName "\"${filter_str}\" ")
filters_snp=$filters_snp$filter0
done
filters_snp=$(echo $filters_snp |perl -pe 's/ $//')

for dat in $(less $tools_path/script/filters_indel.txt); do

filter=$(echo $dat |cut -f2)
filter_str=$(echo $dat|cut -f1)
filter0=$(echo --filterExpression \"$filter\" --filterName "\"${filter_str}\" ")
filters_indel=$filters_indel$filter0
done
filters_indel=$(echo $filters_indel |perl -pe 's/ $//')


echo;echo;echo gatk VariantFiltration
echo java -Xmx$java_memory -jar $gatk -R $ref_genome -T VariantFiltration -o $out_prefix.snp.filtered.vcf --variant $out_prefix.snp.af.vcf $filters_snp |tee .snp.sh

. .snp.sh
	

echo;echo;echo gatk VariantFiltration
echo java -Xmx$java_memory -jar $gatk -R $ref_genome -T VariantFiltration -o $out_prefix.indel.filtered.vcf --variant $out_prefix.indel.af.vcf $filters_indel |tee .indel.sh

. .indel.sh

# echo;echo;echo gatk CombineVariants
# java -Xmx$java_memory -jar $gatk \
	# -R $ref_genome \
	# -T CombineVariants \
	# --variant $out_prefix.snp.filtered.vcf \
	# --variant $out_prefix.indel.filtered.vcf \
	# -o $out_prefix.snp.indel.filtered.vcf \
	# -genotypeMergeOptions UNSORTED
	
	
echo;echo;echo gatk CombineVariants
java -Xmx10g $tmp -jar $gatk \
	-R $ref_genome \
	-T CombineVariants \
	--variant $out_prefix.snp.filtered.vcf \
	--variant $out_prefix.indel.filtered.vcf \
	-o $out_prefix.snp.indel.filtered.vcf \
	-genotypeMergeOptions UNSORTED


echo;echo;echo snpsift annotate id
java -Xmx10g $tmp -jar $snpsift \
annotate \
-id $data_path/ncbi/dbsnp/All_20150605.vcf.gz \
$out_prefix.snp.indel.filtered.vcf \
> $out_prefix.snpsift.dbsnp.vcf


echo;echo;echo snpsift annotate clinvar
java $tmp -jar $snpsift \
annotate \
-id ${data_path}/ncbi/clinvar/clinvar_20150106.vcf \
$out_prefix.snpsift.dbsnp.vcf \
> $out_prefix.dbsnp.clinvar.vcf




echo;echo;echo snpsift annotate cosmic
java $tmp -jar $snpsift \
annotate \
${data_path}/cosmic/CosmicCodingMuts.vcf \
$out_prefix.dbsnp.clinvar.vcf \
> $out_prefix.id.cosmic.vcf


echo;echo;echo snpsift annotate varType
java $tmp -jar $snpsift \
varType \
$out_prefix.id.cosmic.vcf \
> $out_prefix.cosmic.var_type.snpsift.vcf

echo;echo;echo gatk VariantAnnotator
java $tmp -jar $gatk \
	-R $ref_genome \
	-T VariantAnnotator \
	--variant $out_prefix.cosmic.var_type.snpsift.vcf \
	-o $out_prefix.cosmic.var_type.vcf \
	-L $interval \
	-A VariantType


echo;echo;echo snpsift annotate gwasCat
java $tmp -jar $snpsift \
gwasCat \
-db ${data_path}/snpeff/gwascatalog.txt \
$out_prefix.cosmic.var_type.vcf \
> $out_prefix.var_type.gwas.vcf


echo;echo;echo snpsift dbnsfp
java $tmp -jar $snpsift \
dbnsfp \
-db ${data_path}/dbnsfp/dbNSFP2.9.txt.gz \
-f Ensembl_transcriptid,Uniprot_acc,Interpro_domain,SIFT_score,Polyphen2_HVAR_score \
$out_prefix.var_type.gwas.vcf \
> $out_prefix.gwas.dbnsfp.vcf

# mv $out_prefix.var_type.gwas.vcf $out_prefix.gwas.dbnsfp.vcf
	
echo;echo;echo snpeff ann snpeff2gatk
java $tmp -jar $snpeff \
	-c ${snpeff_path}/snpEff.config \
	-v \
	-lof \
	-noStats \
	$snpeff_db_version \
	$out_prefix.gwas.dbnsfp.vcf \
	> $out_prefix.dbnsfp.snpeff2gatk.vcf
	
echo;echo;echo VEP annotation no pick '(can not run in background)'
$tools_path/ensembl-tools-release-82/scripts/variant_effect_predictor/variant_effect_predictor.pl \
--offline \
--force_overwrite \
--dir $data_path/ensembl/ \
--fasta $ref_genome \
--refseq \
--hgvs \
--format vcf \
--vcf \
--input_file $out_prefix.dbnsfp.snpeff2gatk.vcf \
--output_file $out_prefix.vep.vcf




:<<!
echo
echo
echo $gatk VariantFiltration
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o $out_prefix.gatk.filter.vcf \
	--variant $out_prefix.dbnsfp.snpeff2gatk.vcf \
	--filterExpression "QD < 2.0" \
	--filterName "gatk_QD_less_than_2" \
	--filterExpression "MQRankSum < -12.5" \
	--filterName "gatk_MQRankSum_less_than_-12.5" \
	--filterExpression "ReadPosRankSum < -8.0" \
	--filterName "gatk_ReadPosRankSum_less_than_-8" \
	--filterExpression "ReadPosRankSum < -20.0" \
	--filterName "gatk_indel_ReadPosRankSum_less_than_-20" \
	--filterExpression "InbreedingCoeff < -0.8" \
	--filterName "gatk_indel_InbreedingCoeff_less_than_-0.8"
	
!


echo;echo; echo vcf header modify
cat $out_prefix.vep.vcf \
| awk -F', ' '{if (/##INFO=<ID=\w+, /) { print $$out_prefix","$2","$3","$4}else{print $0}}' \
> $out_prefix.header_modified.vcf

# $bgzip -c $out_prefix.header_modified.vcf > $out_prefix.header_modified.vcf.gz
# $tabix $out_prefix.header_modified.vcf.gz

if test "$somatic" = "T"; then
echo filter_high_low.pl
filter_high_low.pl $out_prefix.header_modified.vcf $out_prefix.high.txt $out_prefix.low.txt $out_prefix.header.txt
else
echo filter_pass.pl
filter_pass.pl $out_prefix.header_modified.vcf $out_prefix.high.txt $out_prefix.low.txt $out_prefix.header.txt
fi

cat $out_prefix.header.txt $out_prefix.high.txt > $out_prefix.high.vcf
cat $out_prefix.header.txt $out_prefix.low.txt > $out_prefix.low.vcf

anno_bt.sh $out_prefix.high.vcf
txt2xls.pl $out_prefix.*multianno.txt $out_prefix.anno.high_confidence.xls
format_jxl.sh $out_prefix.anno.high_confidence.xls $out_prefix.anno.high_confidence.jxl.xls

anno_bt.sh -p$out_prefix.2 $out_prefix.low.vcf
txt2xls.pl $out_prefix.2.*multianno.txt $out_prefix.anno.low_confidence.xls
format_jxl.sh $out_prefix.anno.low_confidence.xls $out_prefix.anno.low_confidence.jxl.xls


if test "$single_cell" = "T"; then
	# nprot filter
	nprot_filter.sh $out_prefix.high.vcf

	# impact anno
	impact_anno.sh $out_prefix.nprot.high.vcf
	impact_anno.sh $out_prefix.nprot.low.vcf
else
	impact_anno.sh $out_prefix.high.vcf
fi
impact_anno.sh $out_prefix.low.vcf


. $cmd_done

:<<eof2
echo;echo;echo snpsift filter
java -Xmx$java_memory -jar $snpsift \
filter \
"(FILTER = 'PASS')" \
$out_prefix.header_modified.vcf \
> $out_prefix.high.vcf

# Allele_frequency_less_than_20%;HaplotypeScore_more_than_200 is just a string for snpsift filter X = 'Y'

$bgzip -c $out_prefix.high.vcf > $out_prefix.high.vcf.gz
$tabix $out_prefix.high.vcf.gz

eof2

:<<eof1
echo;echo;echo bcftools view for filter FILTER
$bcftools view \
-f "FP_VARIANTCALLER,FP_PRIMERDIMER" \
$out_prefix.header_modified.vcf.gz \
| bgzip -c \
> $out_prefix.fp.vcf.gz

$tabix $out_prefix.fp.vcf.gz
gunzip -c $out_prefix.fp.vcf.gz > $out_prefix.fp.vcf


$bcftools concat -aD $out_prefix.high.vcf.gz $out_prefix.fp.vcf.gz |bgzip -c > $out_prefix.high_fp.vcf.gz
$tabix $out_prefix.high_fp.vcf.gz

if [ ! -d isec ]; then mkdir isec; fi
$bcftools isec -p ./isec $out_prefix.header_modified.vcf.gz $out_prefix.high_fp.vcf.gz

cp ./isec/0000.vcf $out_prefix.low.vcf
eof1



